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ABSTRACT 

The Method of Light Curve Simulations is a tool that has been applied to X-ray monitoring ob- 
servations of Active Galactic Nuclei (AGN) for the characterization of the Power Density Spectrum 
(PDS) of temporal variability and measurement of associated break frequencies (which appear to be 
an important diagnostic for the mass of the black hole in these systems as well as their accretion state). 
It relies on a model for the PDS that is fit to the observed data. The determination of confidence 
regions on the fitted model parameters is of particular importance, and we show how the Neyman 
construction based on distributions of estimates may be implemented in the context of light curve 
simulations. We believe that this procedure offers advantages over the method used in earlier reports 
on PDS model fits, not least with respect to the correspondence between the size of the confidence 
region and the precision with which the data constrain the values of the model parameters. We plan 
to apply the new procedure to existing RXTE and XMM observations of Seyfert I galaxies as well as 
RXTE observations of the Seyfert II galaxy NGC 4945. 

Subject headings: methods: data analysis — methods: statistical — X-rays: galaxies 



1. INTRODUCTION 

The study of the temporal variability of the 
X-ray flux from accreting black holes has re- 
vealed a com plex behavior of the acc retion 
flow (e.g. iMushotzkv. Done fe Pounds! 119931 : 
iRemillard fc McClintockl 120061) . " One widely ap- 
plied tool for characterizing the variability is the Power 
Density Spectrum (PDS). The shape of the broad-band 
PDS as well as the location of identifiable features such 
as breaks and quasi-periodic oscillations provide the ob- 
servational constraints for physical models of the system 
that generates the variability. Of particular interest is 
the apparent linear scaling between the high-frequency 
break timescale and the black ho le mass in these system s 
over many orders of mag nitude (jMcHardv et al.ll2004h . 

The analysis of the PDS of Active Galactic Nu- 
clei (AGN) is complicated by the question of how to 
assign uncertainties to the Fourier amplitudes mea- 
sured from one observed realiza tion of what is pre- 
sumed to be a stochastic process ([Lawrence et al.lfl987l : 
iMcHardv fc Czernvl H98l . A measure of the expected 
spread in the observed values is essential for the correct 
interpretation of the Fourier spectrum. In addition, in- 
trinsic properties of the Fourier transform (red noise leak 
and aliasing), uneven sampling of the time series, and 
measurement uncertainty in the count rate distort the 
spectrum; these effects need to be corrected for when 
quantifying the shape of the broad-band spectrum. A 
method based on Monte Carlo simulations to determine a 
reliable measure of the PDS uncertainties an d to account 
for th ese distortions was first proposed by iDone et all 
(1992). The main feature of the method is the concept 
of simulating the possible range of realizations of the un- 
derlying process and incorporating the shape-distorting 
effects of uneven sampling, red noise leak, and aliasing. 

1 Kavli Institute for Particle Astrophysics and Cosmology, 
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Uncertainties on the Fourier spectrum are determined 
using light curves generated from a model for the PDS, 
and the level of agreement between the model and the 
data is quantified by a \ 2 fit statistic. By analogy to X- 
ray spectral fitting, the application of these observational 
effects to the chosen model results in the "folded model" , 
which is then used in the comparison to the observations. 

Sub sequent work by lUttlev. McHardv fc Papadakis! 
(2002]) (hereafter UMP02), incorporating the recom- 
mendations for the simul ation of stochastic processes in 
iTimmer fc Koenia (|1995l ) . led to a canonical method for 
the analysis of AGN X-ray light curves, obtained mainly 
from RXTE and XMM-Newton. The process to be mod- 
eled is expressed in the form of a parametric expression 
for the PDS; depending on the complexity of the model, 
a varying number of adjustable parameters determine 
the shape and normalization of the model PDS. In 
addition to the updated Monte Carlo simulations, the 
authors present detailed procedures for the statistical 
evaluation of the model fit, i.e. the assignment of a 
goodness-of-fit measure and the derivation of confidence 
regions on model parameters. The fit statistic, which 
was dubbed the "rejection probability" by subsequent 
authors, is now different from a standard x 2 statistic. 
This development toward a statistically more sophisti- 
cated technique was influenced by considerations about 
the resolution of the PDS, which often needs to be 
compromised to satisfy the con ditions under which the 
X 2 st atistic may be used safely (jPapadakis fc Lawrence! 
Il99l . This new method has found widespread appli- 
cability; the result s reported inj Markowitz e t al.l (l2003f) 



herea fte r M031. IMcHardv et alJ (l2004h. IMarkowitz 



200l) . IMcHardv et all (l2005lh lUttlev fc McHardv 



(2005), IMcHardv et al.l (12007D. iSummons et al. 
(120071). lArevalo. McHardv fc Summons! &20oW. and 
iMarshall. Rvle fc Miller! (|2008!) are all bas ed on it. Our 
initia l report on the PDS of NGC 4945 (|Mueller et all 
l2004h similarly took the published method and in- 
troduced some additional changes. (In contrast to 
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the above papers. iGreen, McHardvfeDon^] (TT9 99) , 



IVaughan. Fabian fc Nandral (120031 ). IVaughan fc Fabianl 
(2003). and lAwaki et all ()2005l ) implement Monte Carlo 
simulations for the derivation of uncertainties on the 
PDS, but use the standard \ 2 statistic for the model 
fit.) 

In the general case of fitting a model to a set of data, 
the derivation of best fit values of the model parame- 
ters (called "point estimation" in statistics) involves the 
identification of the location in parameter space at which 
the fit statistic attains an extremum 2 . These best fit val- 
ues are called "estimates"; the recipe for finding an es- 
timate for a particular parameter is called the "estima- 
tor" . Point estimates by themselves are of limited use. 
Instead, confidence regions on the fitted model parame- 
ters characterize how well the data constrain the model, 
and goodness-of-fit tests may be applied to test whether 
the chosen model is an adequate description of the data 
and whether certain models are favored over others. 

The definition of a confidence region is crucial to its 
proper interpretation. A confidence region (with associ- 
ated confidence level C) is a region in parameter space 
computed from the measured data that has a probabil- 
ity C of containing the true set of parameter values. In 
other words, if the measurement were repeated, a differ- 
ent confidence region would be obtained for each data set, 
but a fraction C of them would enclose the (unknown) 
point in parameter space on which the measured data 
are based. This of course assumes that the model under 
consideration is the correct one; if a goodness-of-fit test 
indicates that the chosen model is a bad description for 
the data, then confidence regions on its parameters are 
of little value. 

Confidence regions are often interpreted as expressing 
the precision with which the model parameters may be 
determined given the data. In practice, for any given fit- 
ting procedure, there are usually several plausible meth- 
ods that produce regions with the required property to 
make them confidence regions. A useful additional con- 
sideration is therefore whether the size of the region that 
a chosen method returns depends appropriately on the 
measurement uncertainties. Furthermore, the value of 
the fit statistic at the location of the best fit should have 
no or only a weak influence on the size. 3 

We review some of the concepts of model fitting using 
the x 2 statistic in Section [2] and we show how, in a sim- 
ple toy model set-up, the A% 2 prescription for finding 
confidence regions satisfies the consideration above. The 
"rejection probability" as a fit statistic is evaluated on 
the same criteria in Section[3J and the strong dependence 
of the size of the confidence region on the minimum re- 

2 The use of Monte Carlo simulations for the derivation of the 
folded model in the case of PDS fitting results in a fit statistic that 
can not be expressed in closed form as a function of the parameters. 
Numerical methods therefore need to be employed to search for the 
location of the extremum. 

3 By way of example, in a standard x 2 fit, under the assumption 
that the chosen model is the correct one and that there are no 
systematic errors in the measurement, the minimum value of x 2 in 
a given fit is fully determined by the ratio of the actual amount of 
statistical fluctuations in the data to the expected amount and does 
therefore not depend on the size of the measurement uncertainties. 
In situations encountered in practice, the conditions under which 
this is true are often violated to a certain degree, such that a small 
influence on the size of the confidence region cannot be ruled out. 



jection probability demonstrated. We then introduce the 
Neyman construction based on simulated distributions 
of estimates in Section 2] as an alternative to the use of 
the rejection probability. This paper does not present 
any actual results obtained from the proposed method. 
However, we outline in Section possible changes that 
may occur if PDS model fits obtained using contours of 
constant rejection probability, including our own work 
on NGC 4945, were re-examined using the Neyman con- 
struction. Section [5] summarizes the paper. 

2. POINT ESTIMATION AND CONFIDENCE REGIONS 
USING THE x 2 STATISTIC 

The most familiar fit statistic in Astrophysics is with- 
out a doubt the \ 2 statistic. It applies well to problems 
where the measured quantities are Gaussian distributed 
with known uncertainties. Even in cases where that con- 
dition is not satisfied, the \ 2 statistic can sometimes still 
yield useful parameter estimates. However, its main at- 
traction lies in the ease with which confidence regions 
on fitted parameters can be derived if the distributions 
are Gaussians, namely through the concept of Ay 2 (see 
e.g.lLampton. Margon fe Bowver|[l976HPress et al.lll992l : 
iBevington fe Robinsonll2003h . Any desired significance 
level < a < 1 maps onto a value of A% 2 such that, after 
determining the best fit values of the model parameters 
by minimizing \ 2 1 the region in parameter space bounded 
by the surface of constant A% 2 contains the true set of 
parameter values with a confidence C = 1 — a. 

Let us illustrate the A% 2 procedure on a toy model 
setup to introduce additional notation that we will refer 
back to in subsequent sections. 

Let y be a physical variable that is expected to be 
proportional to a single independent variable x. As part 
of an experiment, y is measured for a fixed set of non- 
equal Xi (i — 1,...,N). The measurement is expected 
to result in Gaussian uncertainties on y, with a constant 
standard deviation a independent of i. Let {yi} (i — 
1, N) be the set of measurements at the corresponding 
values Xi. 

We now wish to fit these data with a model y = kx. 
The x 2 fit statistic is then a function of the one model 
parameter k and the set of observed values {yi} (all sums 
are over i from 1 to N): 



\k,{y*}) = J2 



{yi -kxif 



(1) 



Minimizing x 2 (fc, {yi}) with respect to k yields the es- 
timate k({yi}) and the minimum value of the fit statistic 

Xmin({yJ) : 



and 



KM) 



{yi -kxif 



(2) 



(3) 



Using these two equations, the expression for the 
change in \ 2 as k is varied evaluates to 
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= ^(fc-fc({ 2/l })) 2 . (4) 

To derive the 68% confidence region (i.e. the "lcr" 
uncertainties) 4 around k{{vi}) (significance a — 0.32), 
we set Ax 2 {k, {yi}) — 1.00. The resulting region satisfies 

\k-k({ yi })\<-^=. (5) 

Since we assumed that the measured yi are in fact well- 
described by the model, then each yi has to be drawn 
from a Gaussian distribution around the true value, i.e. 

Vi ~ g(kti Ue Xi,a), (6) 

where g(a, b) is a Gaussian distribution with average a 
and standard deviation b, and ~ denotes "drawn from." 
The estimate k({yi}) is then drawn from the following 
probability density function: 

k({Vi}) ~ 9 ^fctruc, yg^g J • ( 7 ) 

In other words, if the act of measuring the set {yi} is 

repeated many times, k({yi}) will differ from fc truc by less 

than <t/\/J2 x i 68% of the time. It should now be im- 
mediately obvious that the size of the confidence region 
(Equation [S]) is such that in precisely those 68% of cases 
the confidence region includes k tl - uc , confirming what the 
confidence region was designed to express about the ex- 
periment. 

We have thus confirmed that in this simple setup, the 
Ax 2 prescription produces intervals for the model pa- 
rameter k that satisfy the requirements of confidence in- 
tervals. Furthermore, it can be seen from Equation [5] 
that the size of the confidence interval is proportional 
to the measurement uncertainty a and independent of 
Y^ahn We defer to existing pub lications (specifically 
lLampton. Margon fc Bowverlll976l ) for the extension of 
these results to higher-dimensional parameter spaces. 
We only wish to note here that the independence of 
the size of the confidence region on Xmin i s guaranteed 
through the independence of the distribution of Ay 2 on 
X^nhn, as demonstrated in lLampton. Margon fc Bowverl 
lEm, Appendix IV). 

3. THE REJECTION PROBABILITY 

The measurement of the level of agreement between 
the model and the observed data in the method of light 
curve simulations in UMP02 relies on a statistic called 
by subsequent authors the "rejection probability." It is 
defined analogous to a p- value, with the rejection proba- 
bility being one minus the p-value of the measured Xdist 
fit statistic. Differences in best-fit rejection probability 
between different models are used to favor one model over 

4 The term "lcr" is sometimes used to indicate the 68% confi- 
dence region even if the fit statistic is not \ 2 - I n such cases, the 
standard deviation of the distribution of the parameter may not 
have the same interpretation, but the confidence regions parame- 
terized by the confidence C are always well-defined. 



the others (e.g. a broken power law model compared to 
an unbroken power law model), and regions in parame- 
ter space where the rejection probability is less than a 
certain value (e.g. 90%) are then taken as the confidence 
regions for the fitted model parameters. 

3.1. Confidence Regions from Rejection Probability 

By analogy to x 2 fitting, the UMP02 procedure for de- 
termining confidence regions is equivalent to identifying 
the region in parameter space where x 2 is less than some 
critical value. For the c% confidence region, this critical 
value is simply the c-th percentile of the x 2 distribution 
with the appropriate number of degrees of freedom and 
is thus independent of the minimum value of x 2 obtained 
in the fit. The reason why the authors do not rely on 
percentiles of x 2 distributions for the determination of 
confidence regions is that the effective number of degrees 
of freedom varies with position in parameter space. De- 
ciding on the basis of p- values whether a certain point in 
parameter space is included in the confidence region is 
therefore more robust. 

The region produced in this manner do have the re- 
quired property to make them confidence regions, i.e., 
they include the true value of the parameters with the 
desired probability. However, their sizes depend strongly 
on the value of the fit statistic at the location of the best 
fit. If the minimum rejection probability in a fit is just 
below 90%, the contours of 90% rejection probability will 
be found fairly close around the best fit, leading to the 
erroneous conclusion that the data lend themselves to 
the placement of very precise limits on the model pa- 
rameters. Note that a minimum rejection probability of 
90% does not by itself indicate a bad fit, since there is 
still a 10% chance of obtaining a fit as bad or worse due 
simply to statistical fluctuations; we are therefore justi- 
fied in searching for the confidence region associated with 
the parameters of such a fit. Conversely, if the minimum 
rejection probability is very low, the 90% contours will 
enclose a large area. Furthermore, if the minimum re- 
jection probability is above 90%, there will be no 90% 
confidence region at all. 

We illustrate the inverse correlation between the min- 
imum rejection probability and the size of the resulting 
confidence region schematically in Figure [T] This behav- 
ior is apparent in some of the published results (U MP02, 
M03, iMcHardv et~aT1l2007L ISummons et~aT]|2007h . 

3.2. The Empirical Distribution of x^ ist 

The calculation of the rejection probability relies on an 
approximate determination of the empirical distribution 
of Xdist ■ Since the effective number of degrees of freedom 
is a function of the model parameters, the distribution 
is rightfully calculated separately for each grid point in 
parameter space. However, the Xdist values for the simu- 
lated Fourier spectra are only calculated at their original 
grid point and are not the best fit values found by per- 
forming th e point estimati on on each simulated spectrum 
(see e.g. lPress et al.lfl992l Section 15.6). 

The approximation was most likely introduced due to 
considerations about computing time, because the min- 
imization of Xdist over the parameter space for the hun- 
dreds of simulated spectra that are typically involved in- 
curs a significant computational load. However, given the 
advances in computer power since the original method 
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Fig. 1. — Schematic plot illustrating the dependence of the size of 
the confidence region on the minimum value of the rejection prob- 
ability. The solid and dashed lines are stylized representations of 
the behavior of the rejection probability as a function of a model 
parameter k for fits to two different data sets. Both fits return 
the same estimate for fc, but due to statistical fluctuations, the fit 
indicated by the solid line results in a significantly lower minimum 
rejection probability. The procedure for determining the 68% con- 
fidence interval on k for each fit is indicated by the dotted lines, 
showing the projection of the intersection points between the line 
of 68% rejection probability and the respective parabola onto the k 
axis. Even though this is only a schematic representation, and the 
detailed behavior of the rejection probability as a function of any 
of the parameters in a real fit may be more complicated, the neg- 
ative correlation between the size of the confidence region and the 
minimum value of the rejection probability is expected in general. 

was introduced, the reason for the approximation may 
well have fallen away. 

4. CONFIDENCE REGIONS AND GOODNESS-OF-FIT TEST 
USING SIMULATED DISTRIBUTIONS OF ESTIMATES 

We propose the following set of procedures, most im- 
portantly the Neyman construction based on distribu- 
tions of estimates for the derivation of confidence regions, 
as an alternative to the use of the rejection probability 
for PDS model fits. The new method returns confidence 
intervals whose size has the desired property of being in- 
dependent of the value of the fit statistic at the location 
of the best fit. Furthermore, it deals very naturally with 
biased estimators 5 . 

Throughout this section, it is assumed that the Xdist n ^ 
statistic can be calculated for an arbitrary point in pa- 
rameter space through the use of simulated light curves. 
The procedures are however not specific to the Xdist n * 
statistic; any other statistic which attains an extremum 
at the location of the best fit may be substituted for Xdist • 

4.1. Point Estimation 
Xdist mav b e used directly for point estimation, i.e. the 
estimates o bs for the parameters of the model used to 
describe the observed Fourier spectrum Pobs^) are the 
values of the parameters at the grid point that minimize 
Xdist ■ The estimates 6 s i m for any of the simulated light 
curves (used further below) can be found similarly by 

5 Biased estimators are estimators with an expectation value 
different from the true one. The k estimator used in Section [2] is 
unbiased because its expectation value is ktmc (Equation [7} . In 
more complicated situations, such as the PDS fits under considera- 
tion here, one does not generally know a priori whether the chosen 
estimators are biased or not. 



substituting the simulated spectrum in place of the ob- 
served spectrum and minimizing Xdist over the parameter 
space. 

4.2. Goodness-of-Fit and Hypothesis Testing 

In order to test whether the minimum Xdist value of the 
observed Fourier spectrum signifies an acceptable fit, we 
use the simulations to determine the distribution from 
which Xdist i s drawn. The null hypothesis is that the 
measured Fourier spectrum was in fact produced by the 
model under consideration. Let 6b G st be our best guess 
for the true values of the parameters, i.e. the grid point 
closest to the center of the confidence region. (If the es- 
timators are unbiased, Obest can be set equal to 9.) For 
each of the simulated light curves generated for Obest, we 
record its best fit Xdist (already found above in the deter- 
mination of the confidence region) . The goodness-of-fit 
measure is then the familiar p- value of the observed spec- 
trum's minimum Xdist compared against this distribution 
of simulated Xdist values. As such, it expresses the prob- 
ability that a Xdist vahie at least as high as the measured 
one would be obtained by chance; a p- value smaller than 
the desired significance level (e.g. 5%) indicates that the 
null hypothesis can be rejected. 

The reason for choosing Gbost over any other grid point 
is that the distribution of the minimum values of Xdist 
may depend on O. In standard x 2 fitting, the distribu- 
tion of Xmin is independent of 6, being in fact the x 2 
distribution with the appropriate number of degrees of 
freedom. In the framework of Fourier spectral fits, this 
independence appears to be broken, such that the ef- 
fective number of degrees of freedom is a function of the 
model parameters, plausibly because the degree to which 
adjacent bins in the Fourier spectrum are correlated de- 
pends on the amount of red noise leak (Mueller et al. , in 
preparation). Using Gbost ensures that the Xdist distribu- 
tion thus found approximates as closely as possible the 
one from which the measured Xdist was m ^ ac ^ drawn. 

Note that, up to the approximation to the distribution 
of Xj ist used by UMP02, this procedure is essentially 
equivalent to the calculation of the rejection probabil- 
ity. The p- value is however only used as a goodness-of-fit 
measure and not for finding the confidence intervals. 

If more than one model are under consideration to ex- 
plain the measured data, e.g. when one would like to test 
for the presence of a break in the Fourier spectrum, a de- 
cision statistic for hypothesis testing needs to be set up. 
In the framework of the Xdist n * statistic, the difference 
in best fit Xdist values between two models is a natural 
choice for such a statistic (by analogy to the F-test for 
the x 2 fit statistic). The simulations can once again be 
used to determine the distribution of this difference, from 
which the critical value corresponding to a desired power 
of the test ( "statistical significance" ) may be derived. We 
do not further elaborate on this procedure here, since the 
numbers and decisions involved depend on a balance be- 
tween the sensitivity and specificity of the test that can 
only be calculated using actual simulations. 

4.3. Confidence Regions 

We implement the Neyman construction (|Neymanl 
Il937f ) based on simulated distributions of estimates to 
find confidence limits on model parameters: Let C be the 
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Fig. 2. — Schematic illustration of the Neyman construction applied to a model fit with two adjustable parameters 8± and 82- The plot 
on the left shows elliptical regions obtained from the distribution of the estimates that encompass an unspecified, but constant, fraction 
C of the distribution. The corresponding trial values of the parameters for each ellipse are indicated by the crosses. The offset between 
the crosses and the centers of the ellipses implies a bias in the estimators, kept constant as a function of the parameters in this simple 
example. The location of the observed best-fit (61, O2) is denoted by the cross-hairs. The solid ellipses include the observed best-fit values 
of the parameters, the dashed ones do not. By the prescription of the Neyman construction, the parameter values associated with the solid 
ellipses are added to the confidence region, the others are not. The plot on the right shows the elliptical confidence region (confidence = 
C) that would be obtained if this procedure were to be repeated for all possible trial values of the parameters. The observed best-fit values 
are once more indicated by the cross-hairs. Note how the estimator bias identified earlier results in a confidence region whose center is 
offset from the observed best-fit values. 



desired confidence, e.g. 68% or 90%, and o b s the esti- 
mates for the observed Fourier spectrum as found above. 
Consider now an arbitrary grid point in parameter space, 
©trial- Using the simulated light curves generated for 
that point, we can determine the distribution of esti- 
mates ©sim and derive a region in parameter space that 
encloses a fraction C of them. If © bs is inside that re- 
gion, ©trial is included in the confidence region, otherwise 
it is not. Figure [5] shows this graphically in an imagined 
two-parameter model fit. 

It is easy to see that the size of the confidence regions 
thus obtained is independent of the value of the minimum 
Xdisti n0 information about the measured data enters the 
calculation of the distribution of estimates, and only the 
estimates o bs are used in the subsequent mapping of the 
confidence region. Also, the distribution of estimates, be- 
ing a measure of the possible range of parameter values 
that might be observed given the design of the measure- 
ment, will be broad or narrow depending on the uncer- 
tainties in the observed Fourier spectrum. Therefore, 
the size of the confidence region will scale with the un- 
certainties, preserving the intended interpretation that 
the size of the confidence region expresses the precision 
with which fitted model parameters can be determined. 

As a consequence of using distributions of estimates 
in the Neyman construction, biased estimators will be 
corrected, and the region returned by the algorithm has 
the required property of enclosing the true values of the 
parameters with probability C. In principle, if any of 
the estimators are strongly biased, © bs might lie outside 
the confidence region. The confidence region is however 
always more meaningful in the determination of probable 
values of model parameters than a (possibly biased) point 
estimate. 

Note also that the shape of the region enclosing the 
fraction C of all simulations may be freely chosen by the 
user. This freedom is intrinsic to the Neyman construc- 
tion. However, in order to obtain the tightest possible 
constraints on the parameters, the smallest region should 



be used. 

A complication arises out of the limitation of being able 
to evaluate Xdist 0Iu y on a grid in parameter space. The 
distributions of the estimates are therefore composed of 
finite-size volume elements centered on each grid point. 
The use of smoothed approximations to the empirically 
derived distributions can eliminate this step-wise behav- 
ior; for a one-dimensional parameter space, if the dis- 
tribution of the estimate is sharply peaked, one might 
use a Gaussian fit and determine from the fitted average 
and standard deviation whether the trial point lies inside 
the confidence interval. Similarly, in two dimensions, the 
use of ellipses fitted to the empirical distribution to en- 
close the required fraction of simulated estimates might 
be appropriate. 

Let us illustrate the procedure of finding confidence 
regions using simulated distributions of estimates on the 
toy model introduced in Section [2j The difference to 
the procedure described above is that the estimate fc can 
be found analytically and does not rely on a grid search 
using simulations. However, let us suppose that sim- 
ulations were set up for this simple problem. We use 
Equation [5] to calculate fc bs for the measured data {yi}. 
Let fctriai be an arbitrarily chosen real number, for which 
we want to determine whether it is inside the confidence 
interval. Using simulations with fctriai as the "true" pa- 
rameter value (i.e. randomizations of the observations 
as given by Equation [6|) , we would then find (Equation 
EJ) that the probability density function of the estimates 
for this trial value is given by g (fecial , a/ x i ) ; i- e - 
a Gaussian distribution centered on fctriai- (In this toy 
model, the distribution of the estimates is translationally 
invariant under changes in fctriai j this feature is not ex- 
pected in general.) The smallest interval enclosing 68% 
of the simulations is comprised of the values within one 
standard deviation from fctriai, and by the prescription of 
the method fctriai is included in the confidence region if 
and only if 
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Atrial - fcobsl < /=S=S ) ( 8 ) 

thus recovering the confidence region in Equation [5j 

4.3.1. The Confidence Interval for the Model Normalization 

The model to be fit to the data usually includes an 
overall normalization factor that carries through to the 
model prediction as a multiplicative factor. In this situ- 
ation, the derivation of confidence intervals on the model 
normalization can be simplified. In practice, simula- 
tions need only be done once for an arbitrary normaliza- 
tion, since the model prediction and uncertainties for any 
other normalization may be calculated simply by scal- 
ing. (For a discussion of the complications introduced 
by measurement uncertainties, see Appendix IA.1I ) For 
point estimation, the best fit normalization at any point 
in parameter space 6 can easily be found, either analyt- 
ically or through a numerical search. This procedure is 
unchanged from UMP02 (section 4.2). The estimates for 
the remaining parameters, either for the measured data 
or for any of the simulated spectra that may be substi- 
tuted for it, are then once again the values of the param- 
eters at the grid point in the remaining parameter space 
where Xdist attains a minimum. 

For the derivation of confidence regions, we again wish 
to determine whether a certain grid point in parameter 
space, with normalization Atrial an d remaining parame- 
ters ©trial, is included at a given significance level. Let 
No be the original normalization with which the light 
curves at tr iai were simulated, and /o(-A, @) the cor- 
responding probability density function of the estimate 
distribution with its dependence on the normalization 
estimate N and the estimates of the remaining model 
parameters 0. Because of the multiplicative nature of 
the model normalization, this distribution becomes scale- 
invariant along the N axis (see Appendix [AJ , such that 
if a different normalization had been used to simulate 
those light curves, the distribution would simply be an 
appropriately scaled version of fo(N, 0). 

The point estimation on the observed Fourier spectrum 
defines the location of the best fit, given by (N b s , ©obs)- 
We now determine the region R in the (N, 0) space that 
encloses the required fraction C of the estimate distri- 
bution (e.g. 68%). Now consider the line in parameter 

space along which = © ob s, i.e. the line parallel to 
the N axis that passes through the best fit. This line 
either intersects the boundary of R in a finite number of 
points, or else no intersection points exist. The region R 
is in most cases convex, such that there are either zero or 
two intersection points. We will only consider these two 
cases here; the procedure is easily generalized to non- 
convex regions that may result in additional intersection 
points. 

In the case of zero intersection points, (Atrial, ©trial) 
is excluded from the confidence region for all values of 
Atrial- In the other case, let us denote the values of the 
normalization at the two intersection points by A^jow 

6 The parameter space now under consideration excludes the 
model normalization as a parameter, due to its special treatment. 



and A^high- Because of the scale-invariance of the esti- 
mate distribution, these values are proportional to the 
original normalization A^o, such that, for any other nor- 
malization A" that could have been used to generate the 
light curves, 

N ~ 

N ow (N) = — Ao,iow (9) 

and 

N » 

N high (N) = — N 0Mgh . (10) 
i\ 

The condition for (Atrial, ©trial) to be included in the 
confidence region now reduces to whether the observed 
value of the normalization is located between these two 
bounds, i.e. 

Mow (Atrial) < A> obs < A> high (AT tria i), (11) 

which is equivalent to the condition on Atrial 

AT -^<iV trial <iVo-^. (12) 

Ao,high A'ojow 

In summary, the estimate distribution found from light 
curves simulated at ©trial, with normalization N , may 
be used to derive the bounds Ao,iow and Ar^high, from 
which the limits on Atrial (for a given ©trial) may be 
calculated. The confidence region in the full (N, ©) pa- 
rameter space finally may be mapped out by repeating 
the procedure for different values of ©trial- 

5. DISCUSSION 

Applying the criterion whether the procedure for de- 
termining confidence regions returns regions whose size 
scales appropriately with the uncertainties in the data, 
we believe that the Neyman construction based on simu- 
lated distributions of estimates offers a viable and advan- 
tageous alternative to the use of the rejection probability. 
While neither UMP02 nor the authors of the subsequent 
papers utilizing the method make specific claims regard- 
ing the statistical properties of the regions obtained, any- 
one not familiar with the data analysis at a sufficient level 
of detail will tend to interpret the quoted uncertainties 
on the best-fit values of the model parameters as indica- 
tive of the precision with which the data constrain those 
values. 

We wish to stress however that this does not imply 
that previously reported results are inherently flawed. 
The confidence limits on break frequencies and power 
law indices for AGN PDS fits may turn out to be dif- 
ferent under the application of the new method, but it 
remains to be seen whether any of these changes are large 
enough to substantially change the interpretation of the 
observations. Specifically, we do not expect that the lin- 
ear scaling between br eak timescale and black hole mass 
(jMcHardv et al.ll2004D will be affected even if the con- 
fidence intervals on some of the data points were to be 
modified. 

It is likely that the precision with which the break fre- 
quencies for Fairall 9, NCG 4151 (both from M03), and 
MCG-6-30-15 (UMP02) have been reported is too opti- 
mistic. Similarly, the confidence regions for the peak fre- 
quencies of the Lorentzians in the fit to the PDS of Ark 



7 



564 may be too small, especially considering that even a 
small increase in the size of the confidence contour in a 
plot where the axes are the logarithms of the peak fre- 
quencies (Figure 9 in their paper) has a disproportionate 
effect on the uncertainties on the frequencies themselves. 
On the other hand, some break frequencies may in fact 
be better determined with current data than reported 
in the literature. Examples of fits where the minimum 
rejection probability turned out to be particularly low 
include NGC 5548 and NGC 3516 (M03). 

A related issue is the use of contours of constant 
rejection probability as confidence limits on combina- 
tions of paramete rs, such as the ratio o f break frequen- 
cies (Figure 11 inlMcHardv et al.ll2005l and Figure 6 in 
lUttlev & McHardvl 120051) . " Limits on this ratio are used 
as key pieces of evidence to motivate the association of 
the AGNs under consideration (NGC 3227 and MCG-6- 
30-15) with the analogue of the high/soft accretion state 
in galactic black hole X-ray binaries. Given that the 
confidence regions were derived using the rejection prob- 
ability, the quoted confidence values with which certain 
ranges of ratios are excluded in those reports may or may 
not in fact be supported by the data. We do not expect 
that the use of the new method will alter the general 
direction of these results, i.e. that the ratio of break 
frequencies in these AGN is likely to be higher than ex- 
pected for the low/hard state, but a re-analysis of the 
observations focusing on the doubly-broken power law 
model might be warranted. The calculation of the sta- 
tistical significance with which the model where the ratio 
of these break frequencies is fixed at a value of 30 may be 
rejected in favor of the original model where both break 
frequencies are allowed to vary forms an additional test 
on these data. If both of these lines of evidence produce 
mutually consistent results, the case for the classification 
of these AGN as analogues of galactic X-ray binaries in 
the high/soft state will be strengthened. 

On the question of the statistical significance of breaks 
in the PDS of AGN, we believe that additional work is 
needed to make the values that have been reported more 
secure. Table 5 in M03 lists the quantity Act that was de- 
signed to express the increase in likelihood of the fit once 
a break is added to the PDS model. It is however not 
clear from the description whether Act was calculated 
using the rejection probability or the underlying Xdist 
values. As outlined in Section l4~2l the Xdist ^ statistic 
does lend itself to the formulation of such a hypothesis 
test. A validation of critical values of differences in Xdist 
and their corresponding statistical significances will be 
required. This includes using the simulated light curves 
to calculate type I and type II errors (rate of false posi- 
tives and false negatives) or, equivalently, the specificity 
and sensitivity of the hypothesis test. As far as we are 
aware, the amount and quality of data needed to effect 
a reliable detection of a break at a significance level of 
5%, say, is an unanswered question. A systematic inves- 
tigation in this area, using both real and simulated data, 
might uncover general considerations that would be in- 
valuable for the design of future observatories for AGN 
timing research. 

The Monte Carlo method for calculating folded mod- 
els to include observational effects may have applications 
outside of PDS fits to AGN X-ray light curves. In par- 



ticular, the associated procedure of finding confidence 
limits on fitted parameters using simulated estimate dis- 
tribution could be used for X-ray or 7-ray spectral fits 
in the low counts-per-bin limit, where the discrete na- 
ture of the Poisson process becomes important. The 
study of the PDS of galactic X-ray binaries may ben- 
efit from an application of the Monte Carlo method as 
well. The shorter time scales for characteristic varia- 
tions and the extensive archive of X-ray observations 
allow for a much more detailed investigation into the 
shape of the broad-band variability spectrum, includ- 
ing the direct observation of independent realizations 
of the underlying proce ss (e.g. iPottschmidt et al.l [20031 
iDone fc Gierlinskll2005f ). The measurement of the dis- 
tribution of power in individual frequency bins is of par- 
ticular interest in this case. Competing physical mod- 
els for the variability in these sources predict distinc- 
tive properties of the stationarity and degree of stochas- 
ticity of the proces s unde rl ying the observation s (e.g . 
i Poutanen fc Fabianl 119991: iMaccarone fe CoppI 120021: 
iMinutti fc Fabian! 120041: lUttlev. McHardv fc VaugharJ 
120051 : IZvcki fc Maciolek-Niedzwieckil |2005[ ). Adopting 
the Monte Carlo simulations for the analysis of galac- 
tic X-ray binaries, specifically the comparison between 
predicted and observed distributions of the Fourier am- 
plitude, may lead to tests of certain elements of these 
models. Furthermore, tools beyond the PDS for the in- 
vestigation of these kinds of stochastic proc esses, such 
as the bispectrum (jVaughan fc Uttlevl I2007T ). are more 
sophisticated in their treatment of the underlying vari- 
ability process, but they will also continue to, at least for 
a while, produce results that are not as validated in their 
interpretation as those derived from standard Fourier 
analysis. Monte Carlo simulations will likely remain es- 
sential for the important comparison of these tools to the 
PDS; a solid statistical foundation is in turn essential for 
these simulations. 

6. CONCLUSION 

Evaluated by the criterion whether the sizes of the con- 
fidence regions express the precision with which the data 
constrain the model parameters, we have shown that the 
use of simulated distributions of estimates (Section I4.3[) 
is preferable to the rejection probability (Section l3.1jl . 
Confidence regions determined from the latter do have 
the required property of enclosing the true values of the 
parameters with the given probability, but their size is 
highly variable depending on the minimum value of the 
fit statistic at the location of the best fit. The method 
based in the former is computationally more intensive, 
but is the only way known to us to derive meaningful 
uncertainties on fitted model parameters in the absence 
of a better-understood fit statistic. 

The end products of the application of the set of pro- 
cedures in Section |4~T1 through 14.31 to the observations of 
an AGN are the best fit values of the parameters for the 
model-dependent description of the PDS, the associated 
confidence limits, and the goodness-of-fit of the model (p- 
value of the minimum Xdist)- Depending on the nature 
of the investigation, more sophisticated statistical tests 
may be employed to test different hypotheses against the 
same data set, or to quantify the observed variations in 
the parameter values between different AGN. 

We are in the process of applying the new method to 



the RXTE observation of the Seyfer t II galaxy NGC 4945 
for which we reported first results in lMueller et al.1 (|2004l ) 
and plan to re-analyze the existing archival RXTE and 
XMM-Newton observations of Seyfert I galaxies with the 
updated procedure. While we don't expect the conclu- 
sions drawn from the analysis of these observations to 
change significantly, this will put the investigation into 
the shape of the PDS in AGN on a statistically more 
solid foundation and make the interpretation of the re- 
sults easier. 
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APPENDIX 

THE SCALE-IN VARIANCE OF THE ESTIMATE DISTRIBUTION FOR THE MODEL NORMALIZATION 

The special treatment of the model normalization in the derivation of confidence regions relics on a property of the 
estimate distribution under the conditions mentioned in the text (Section l4.3.ip . namely that the normalization carries 
through to the model prediction as a multiplicative factor. 

Let Ao be the normalization (hereafter called the "input normalization") that was used to generate a set of light 
curves at an arbitrary point in parameter space Gtriai, and let {P m (vi)} be the Fourier spectrum of one of them, where 
v.i are the frequencies over which the spectrum is measured. Additionally, let (A in , <~> in ) be the estimates for this light 
curve that were found as part of the procedure to determine the estimate distribution (see Section l4~3l for details). 
Because the normalization is an overall multiplicative factor in the generation of these light curves, the {Pi n {^i)} values 
are proportional to No. 

The model to be fit to these data can be written as P(v, N, 9) = N P r (v, 9), where A is the model normalization 
and P r {v, 9) the function describing the dependence of the model on the remaining parameters 9. The folded model 
at the point in parameter space given by N and 9 is summarized in two variables for each frequency bin: the average 
power P s im{vi, A, 9) and the standard deviation AP s i m (^, N, 9) (for details, see Section 4.2 in UMP02). Both of these 
scale with A: 



and 



where P s im,r(^i, 9) and AP s j 



P sim K A, 9) = NP sim , r {ui,e), 

AP sim (v h A, 9) = NAP sim , r (vi, 9), 
• {yi, 9) constitute the folded model for N = 1. The fit statistic 

'i\Jyi) - Np~{vi,®y ° 



xLwe.{W}) = E 



AAP sim , r (^,9) 



(Al) 



(A2) 



(A3) 



is invariant under changes in the input normalization Ao — > r\ Nq (t] > 0) if the same multiplicative factor is applied 
to the model normalization A. As a consequence, since N- m and 9i n are the estimates for this simulated light curve for 
rj = 1, then (n A; n ) and 9j n would have been the estimates if the original normalization had been different by a factor 
77. This applies to all simulated spectra; therefore the distribution of the estimates (Mn, 9j n ) will be scale-invariant 
along the A axis: Let /o(A, 9) be the estimate distribution for the original input normalization Ao (i.e. r\ = 1). For 
any other value of 77, the estimate distribution is then 



me)= 1 -fo f,e 



(A4) 



Note that the above does not require that the normalization be uncorrelated with the other model parameters. The 
invariance of the fit statistic is preserved even if such correlations exist. 

Influence of Measurement Uncertainties 

In the context of PDS model fits, the measurement uncertainties in the light curve manifest themselves as an 
additional noise component in the Fourier spectrum (Poisson level). The scaling of the model prediction with the 
normalization factor is only approximate in this case, since the Poisson level is constant and does not scale with the 
model normalization. However, the intrinsic variability in the light curve by design usually dominates over the Poisson 
level. The confidence interval on the model normalization derived while ignoring this complication is therefore expected 
to approximate closely the more correct one that would be obtained through the usual prescription of simulating light 
curves with different normalizations and deriving the distribution of the estimates in each case. 
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Applicability to PDS with Logarithmically Averaged Power 

In the canonical method of UMP02, -P a i rn(^i;) is actually the average o f the logarithm of the periodogram power, 
which is motivated by the considerations inlPapadakis ^"Lawrence! (|1993f ) . The logarithm of the model normalization 
TV therefore enters the model prediction as an additive constant, while the uncertainties AP s ; m (fj), being standard 
deviations on what are now logarithmic power values whose spread is unaffected by the model normalization, are 
independent of N. The estimate distribution is then translationally invariant along the log N axis, and the expression 
for the bounds on Atrial turns out to be the same as for the linear case (Equation [T2l) . 

Note however that different numerical values for these bounds may be obtained depending on whether the estimate 
distribution is expressed as a function of N or log N . The shape and extent of the region R encompassing the 
desired fraction of the estimate distribution may vary; the smallest such region for example will in general be different 
depending on the choice of variables. In a complete description of the analysis method, it will be important to state 
which variable was used. 
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